home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Business Shareware
/
Business Shareware.iso
/
start
/
hobby
/
aa_51
/
precess.c
< prev
next >
Wrap
C/C++ Source or Header
|
1993-02-09
|
6KB
|
239 lines
/* Precession of the equinox and ecliptic
* from epoch Julian date J to or from J2000.0
*
* Program by Steve Moshier.
*
*
* IAU Coefficients are from:
* J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
* "Expressions for the Precession Quantities Based upon the IAU
* (1976) System of Astronomical Constants," Astronomy and
* Astrophysics 58, 1-16 (1977).
*
* Newer formulas that cover a much longer time span are from:
* J. Laskar, "Secular terms of classical planetary theories
* using the results of general theory," Astronomy and Astrophysics
* 157, 59070 (1986).
*
* See also:
* P. Bretagnon and G. Francou, "Planetary theories in rectangular
* and spherical variables. VSOP87 solutions," Astronomy and
* Astrophysics 202, 309-315 (1988).
*
* Laskar's expansions are said by Bretagnon and Francou
* to have "a precision of about 1" over 10000 years before
* and after J2000.0 in so far as the precession constants p^0_A
* and epsilon^0_A are perfectly known."
*
* Bretagnon and Francou's expansions for the node and inclination
* of the ecliptic were derived from Laskar's data but were truncated
* after the term in T**6. I have recomputed these expansions from
* Laskar's data, retaining powers up to T**10 in the result.
*
* The following table indicates the differences between the result
* of the IAU formula and Laskar's formula using four different test
* vectors, checking at J2000 plus and minus the indicated number
* of years.
*
* Years Arc
* from J2000 Seconds
* ---------- -------
* 0 0
* 100 .006
* 200 .006
* 500 .015
* 1000 .28
* 2000 6.4
* 3000 38.
* 10000 9400.
*/
/* Precession coefficients taken from Laskar's paper: */
static double pAcof[] = {
-8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
-0.235316, 0.07732, 111.1971, 50290.966 };
/* Node and inclination of the earth's orbit computed from
* Laskar's data as done in Bretagnon and Francou's paper:
*/
static double nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10,
-3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
-0.042078604317, 3.052112654975 };
static double inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
-5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
0.002278495537, 0.0 };
extern double J2000; /* = 2451545.0, 2000 January 1.5 */
extern double STR; /* = 4.8481368110953599359e-6 radians per arc second */
extern double coseps, sineps; /* see epsiln.c */
/* Subroutine arguments:
*
* R = rectangular equatorial coordinate vector to be precessed.
* The result is written back into the input vector.
* J = Julian date
* direction =
* Precess from J to J2000: direction = 1
* Precess from J2000 to J: direction = -1
* Note that if you want to precess from J1 to J2, you would
* first go from J1 to J2000, then call the program again
* to go from J2000 to J2.
*/
extern int epsiln();
int precess( R, J, direction )
double R[], J;
int direction;
{
double sinth, costh, sinZ, cosZ, sinz, cosz;
double A, B, T, Z, z, TH, pA, W;
double x[3];
double *p;
int i;
double sin(), cos(), fabs();
if( J == J2000 )
return(0);
/* Each precession angle is specified by a polynomial in
* T = Julian centuries from J2000.0. See AA page B18.
*/
T = (J - J2000)/36525.0;
if( fabs(T) > 2.0 )
goto laskar;
Z = (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR;
z = (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR;
TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR;
sinth = sin(TH);
costh = cos(TH);
sinZ = sin(Z);
cosZ = cos(Z);
sinz = sin(z);
cosz = cos(z);
A = cosZ*costh;
B = sinZ*costh;
if( direction < 0 )
{ /* From J2000.0 to J */
x[0] = (A*cosz - sinZ*sinz)*R[0]
- (B*cosz + cosZ*sinz)*R[1]
- sinth*cosz*R[2];
x[1] = (A*sinz + sinZ*cosz)*R[0]
- (B*sinz - cosZ*cosz)*R[1]
- sinth*sinz*R[2];
x[2] = cosZ*sinth*R[0]
- sinZ*sinth*R[1]
+ costh*R[2];
}
else
{ /* From J to J2000.0 */
x[0] = (A*cosz - sinZ*sinz)*R[0]
+ (A*sinz + sinZ*cosz)*R[1]
+ cosZ*sinth*R[2];
x[1] = -(B*cosz + cosZ*sinz)*R[0]
- (B*sinz - cosZ*cosz)*R[1]
- sinZ*sinth*R[2];
x[2] = -sinth*cosz*R[0]
- sinth*sinz*R[1]
+ costh*R[2];
}
goto done;
laskar:
/* Implementation by elementary rotations using Laskar's expansions.
* First rotate about the x axis from the initial equator
* to the ecliptic. (The input is equatorial.)
*/
if( direction == 1 )
epsiln( J ); /* To J2000 */
else
epsiln( J2000 ); /* From J2000 */
x[0] = R[0];
z = coseps*R[1] + sineps*R[2];
x[2] = -sineps*R[1] + coseps*R[2];
x[1] = z;
/* Precession in longitude
*/
T /= 10.0; /* thousands of years */
p = pAcof;
pA = *p++;
for( i=0; i<9; i++ )
pA = pA * T + *p++;
pA *= STR * T;
/* Node of the moving ecliptic on the J2000 ecliptic.
*/
p = nodecof;
W = *p++;
for( i=0; i<10; i++ )
W = W * T + *p++;
/* Rotate about z axis to the node.
*/
if( direction == 1 )
z = W + pA;
else
z = W;
B = cos(z);
A = sin(z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;
/* Rotate about new x axis by the inclination of the moving
* ecliptic on the J2000 ecliptic.
*/
p = inclcof;
z = *p++;
for( i=0; i<10; i++ )
z = z * T + *p++;
if( direction == 1 )
z = -z;
B = cos(z);
A = sin(z);
z = B * x[1] + A * x[2];
x[2] = -A * x[1] + B * x[2];
x[1] = z;
/* Rotate about new z axis back from the node.
*/
if( direction == 1 )
z = -W;
else
z = -W - pA;
B = cos(z);
A = sin(z);
z = B * x[0] + A * x[1];
x[1] = -A * x[0] + B * x[1];
x[0] = z;
/* Rotate about x axis to final equator.
*/
if( direction == 1 )
epsiln( J2000 );
else
epsiln( J );
z = coseps * x[1] - sineps * x[2];
x[2] = sineps * x[1] + coseps * x[2];
x[1] = z;
done:
for( i=0; i<3; i++ )
R[i] = x[i];
return(0);
}